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LINEAR  STABILITY  OF  SELF-SIMILAR  FLOW: 

8.  IMPLODING  CYLINDRICAL  AND  SPHERICAL  SHOCKS 
IN  THE  C-C-W  APPROXIMATION 


I.  Introduction 

One  of  the  critical  limitations  to  achieving  high  compression  in  a 
spherical  implosion  is  the  degree  of  symmetry  that  can  be  maintained.  This 
in  turn  has  important  implications  for  target  fabrication  techniques  and 
for  laser  or  other  driver  designs .since  it  establishes  the  tolerances 
required  in  the  symmetry  of  these  components. 

An  important  issue  for  understanding  imploding  systems  is  the  stability 
of  a  converging  shock  wave.  This  shock  wave  might  be  used,  for  instance, 
not  only  to  compress  the  fuel,  but  also  to  provide  the  heating  required 
to  create  a  central  ignition  region.  The  final  temperature  achieved  will 
depend  on  how  nearly  spherical  the  shock  wave  remains  during  the  collapse 
process  and  the  shape  of  the  shock  at  the  time  of  reflection. 

A  certain  inherent  stability  of  a  shock  wave  results  from  the  well- 
known  fact  that  a  shock  wave  with  a  smaller  radius  of  curvature  advances 
faster  than  one  with  a  larger  radius  of  curvature.  Thus,  the  part  of  a 
perturbed  shock  front  that  initially  lags  behind  will  accelerate  more 
rapidly  due  to  its  smaller  radius  of  curvature  and  so  will  tend  to  catch  up 
with  the  remainder  of  the  shock  wave.  However,  the  perturbation  may  be  un¬ 
able  to  overtake  the  main  shock,  which  is  accelerating  because  of  convergence, 
or  it  may  be  overdriven,  i.e.,  the  perturbation  may  overshoot  the  stable  position. 

In  order  to  discuss  stability,  it  is  necessary  to  define  what  is  meant 

by  stable  (or  unstable)  behavior.  The  usual  definition  of  stability  in 
Manuscript  submitted  September  16,  1980. 


terms  of  a  growing  or  decaying  mode  amplitude  does  not  adequately  describe 
the  situation  in  imploding  systems.  For  example,  the  amplitude  may  not 
tend  to  zero  as  fast  as  the  average  radius,  or  the  collapse  time  may  be 
of  the  order  of  the  period  of  oscillation  of  the  mode.  A  more  meaningful 
number  for  small-amplitude  perturbations  is  the  rate  of  growth  (or  decay) 
of  the  relative  perturbation  amplitude,  i.e.,  the  ratio  of  the  perturbation 
amplitude  to  the  radius  of  the  zero-order  symmetric  collapse  solution 
(Bernstein  and  Book,  1978). 

Large  initial  perturbation  amplitudes  may  not  decay 

in  the  same  way  as  small-amplitude  perturbations  (Fong  and  Ahlborn,  1979). 

One  can  define  a  radial  instability  in  terms  of  the  maximum  deviation 

of  the  shock  radius  from  the  average  radius  (|R-R  | /R  )•  Another 

1  ave  ave 

kind  of  instability  (kink  instability)  occurs  when  a  small  portion  of 
the  shock  necks  off  from  the  central  region  (Fig.  1).  In  general,  we 
cannot  speak  of  an  imploding  shock  as  being  stable  or  unstable  in  a  clear- 
cut  sense.  Rather,  we  can  ask  whether  it  retains  an  acceptable  degree 
of  symmetry  after  having  collapsed  to  a  volume  sufficiently  small  for 
practical  purposes. 

We  have  developed  an  analytic  and  computational  model  to  investigate 
the  stability  of  converging  shock  waves  in  cylindrical  and  spherical 
geometry.  This  represents  an  extension  to  smaller  radii  of  the  work  of 
Fong  and  Ahlborn  (1979)  on  the  linear  stability  problem.  The  model  equations 
are  described  in  the  next  section,  followed  by  a  linearized  analytic  solution, 
an  account  of  the  numerical  model,  a  comparison  of  the  numerical  and  analytic 
solutions,  an  analysis  of  nonlinear  behavior, and  an  extension  to  problems 
with  a  varying  density  in  front  of  the  shock. 


ig.  1  —  Kink  instabilities  form  when  the  main  shock  accelerates  due  to  spherical  convergence  ahead  of  a 
perturbed  portion  of  the  shock  which  is  unable  to  catch  up  in  spite  of  its  smaller  radius  of  curvature. 


II.  Model  Equations 

The  motion  of  a  converging  shock  wave  can  be  computed  with  great 
accuracy  by  considering  only  changes  in  the  physical  variables  across 
the  shock  front,  ignoring  the  fluid  motion  behind  the  shock  surface.  The 
velocity  of  a  shock  front  in  a  motionless  undisturbed  medium  is  normal  to 
the  front.lt  may  therefore  be  treated  as  a  locally  one-dimensional  motion 
in  a  channel  whose  boundaries  are  determined  by  the  trajectories  of  the 
shock  front.  These  trajectories  form  imaginary  ray  tubes  whose  cross- 
sectional  area  is  related  to  the  Mach  number  by  an  equation  derived  by 
Chester  (1954),  Chisnell  (1955),  and  Whitham  (1957)  (the  C-C-W  approxima¬ 
tion)  .  The  equation  may  be  found  by  substituting  into  the  compatibility 
equation  for  the  characteristic  moving  in  the  direction  of  shock  propagation 
the  fluid  quantities  determined  by  the  Rankine-Hugoniot  relations  across 
the  shock.  Whitham  (1958)  has  shown  that  this  procedure,  the  simplest 
method  of  deriving  the  C-C-W  approximation,  is  equivalent  to  solving  the 
complete  fluid  equations  while  ignoring  the  characteristics  which  travel 
away  from  the  shock  and  are  reflected  back  toward  it.  For  this  reason, 
it  is  most  accurate  for  accelerating  shocks ,where  the  reflected  charac¬ 
teristics  are  unable  to  catch  up  with  the  shock. 

The  result  of  this  model  is  an  equation  for  the  Mach  number  M  of  the 
shock  as  a  function  of  the  cross-sectional  area  A  of  the  ray  tube: 

X(M)MdM  ,  dA  _  A  (1) 

2  +  a  “  u* 

M  -1  A 

vhefe  2 

MM)  -  (2c*l+l/M2)  (l+“  ^=2-),  (2) 

and  a2  -  (T-1)*1  +  2, —  ,  O) 

2oM2  -  (Y-l) 
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where  y  is  the  usual  adiabatic  index.  The  cross-sectional  area  of  the 
ray  tube  may  be  expressed  in  terms  of  the  shock  front  velocity  by  the 
kinematic  relations 

A  “  -A  n  •  (n  x  V)  x  v  (4) 

n  =  -  (I  -  nn)  •  (n  x  v)  x  v;  (5) 

r  =  y  =  n  v,  (6) 

where  (*)  denotes  a  derivative  with  respect  to  time,  n  is  the  unit  vector 
normal  to  the  shock  front  and  r  is  the  surface  location. 

In  numerical  integration,  Eqs.  (1)  and  (5)  are  used  to  propagate 
the  magnitude  and  direction,  respectively,  of  the  shock  front  velocity  in 
terms  of  time,  substituting  A  from  Eq.  (4).  Equation  (6)  then  advances 
the  shock  front  position  by  integration  of  the  velocity.  For  small  per¬ 
turbations  about  symmetric  (cylindrical  or  spherical)  shocks,  however, 
solutions  can  be  obtained  analytically. 
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III.  Linearized  Analytic  Solutions 

Since  all  collapsing  shocks  will  eventually  become  strong,  it  is 
useful  to  analyze  the  strong  shock  limit  M»l.  If  we  assume  that  the 
sound  speed  ahead  of  the  shock  is  constant,  in  addition  to  being  small, 
then  Eqs.  (l)-(3)  reduce  to 

A/A - X  v/v  ,  (7) 

where 


and 


M/a 


1+1  _2yh 

Y  Y-l; 


Using  familiar  vector  identities  we  have 

n  •  (n  x  7)  *  v  ■  -v  7  •  n 
Thus  using  Eqs.  (A)  and  (7)-(10)  we  have 

(X/v)‘  -  V'n-nxCnxV)  •  n  , 

which  only  involves  tangential  derivatives  of  n. 

We  first  look  at  the  zeroth-order  symmetric  solution, 

r  -  v(t); 

n  ■  e  r; 

n  ■  0; 

(X/r)‘  ■  7  •  (r/r)  ■  a/r, 

where  a  -  1  (2)  for  cylindrical  (spherical)  geometry.  Equations 
(12) -(15)  can  be  rewritten  as 

r  (A/r)  -  -A  ^  to  t  .  f  -  ^  to  r“  . 


(8) 

(9) 

(10) 

(11) 


(12) 

(13) 

(14) 

(15) 


Integrating  this  once  yields 


and  finally 


.  ol/\ 

r  r  A  ■  const, 
r-R  (w)X/(i+a> 


(16) 


(17) 

(18) 


where  R  and  u>  are  constants. 


Suppose  we  now  linearize  the  equation  by  assuming  n  =  n^  +  n^,  v  =  v^+  v^, 

r  =  Tq  +  C(rQ,  t) , where  the  subscript  0  refers  to  the  zeroth-order  symmetric 
solution,  and  the  subscript  1  refers  to  a  small  perturbation  about 


this  solution.  From  Eq.(  5  )  we  have 

?1  =  “  (l  '  *(&)>  • 

Expanding  the  gradient  operator 

V  =  VQ  -  (70  f) 
from  Eq.  (11)  we  have  to  first  order 

(-Xvi/vo2)"  =  70  '  ?1 


(t?0  x  V  x  5 


70  +  • 


<v0 


From  Eq.  (19)  after  some  manipulation  we  have 

?i  “  !  •  7o  ?  =  •  (*  -  “oV 


<W 


7o  (?o 


S) 


Here  use  has  been  made  of  the  fact  that  n^  =  e and  the  fact  that 

Vo  =  VW  =  V~oV  /ro 

is  a  symmetric  dyad.  From  Eq.  (6)  we  have 

{  *  "o  V1  +  ~1  V0 

Hence  taking  the  dot  product  with  n^  yields 


5o  ’  ?  =  (?o 

Thus  Eq.  (21)  can  be  written  as 


O' 


v. 


^_|a4)  =  [nr  5  •  V0n0]+  c  •  70  V0  .  nQ 

-  -  7o  *I(i  "  “o5o>*  VV  !)]  +  *  '  7o(-^)- 


(19) 


(20) 


(21) 


(22) 


(23) 


(24) 


(25) 


(26) 


Now  if  we  expand  the  perturbation  in  either  cylindrical  or  spherical 
harmonics  by  assuming  a  sum  of  linearly  independent  terms  of  the  form 


no  ’  C  *  C(t)  cos  (m^)  (cylindrical) 


„m 


S(t)  P„  cos  (mv?)  (spherical)  , 

(6) 


(27) 

(27') 


each  of  which  can  be  treated  independently,  and  define  the  mode-number- 
dependent  coefficient  of  the  Laplacian  by 
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(28) 


Q  *  m 


1(1+1) 


(cylindrical) 
(spherical) , 


(28') 


Eq.  (26)  becomes 


Now  using 


-at/vb‘  =  ?(Q-a)/r 
0  0 


5  = 


we  have 


v  - — f  (X/v  )^~]=ot-  ~J?  r 
0  dr0K  /vO)dt0>  r02  ' 


From  Eq.  (18)  we  have 


d£n  v 


0  -  a 


and  after  a  little  algebraic  manipulation  we  get 

,  d2  ?  ,  a  dC  ,  Q  -  a  . 

A  - - =■  +  -  J  +  C  =  0. 

dr02  ro  dro  ro2 


We  now  seek  a  solution  of  the  form  £  ~  r^  ,  where  6  is  in  general 
a  complex  number  which  satisfies  the  indicial  equation 

X  8(8-1)  +  ctS  +  q-  ot»0,  (34) 


S  -  ^-{X-a  +  [(X-Kt)2  -  4XQ]*5}  .  ( 

Since  we  seek  the  ratio  of  the  amplitude  of  the  disturbance  to  that 

of  the  zeroth  order  solution  we  look  at 

_ .  _  ,8-1  ,  C-(X+a)  +  [(X+a)2  -  4X  ]**}  /2  X  . 

4/^0  rQ  "  rQ  ”  ’ 


Since 


ro  -  c 


X 

X  +  a 


8 


t 


(38) 


C/rQ  - 


-*5  + 


ra-Hx)2  -  4xqi>5/2(x-kx) 


For  the  lowest-order  mode  numbers  8-1  is  real  and  negative, 
indicating  that  the  disturbance  is  geometrically  unstable: 

8  -  1  =  -  ,  0  (39) 

for  2=0  (spherical  coordinates)  or  m=0  (cylindrical  coordinates) ,  and 

6  -  1  *  -  ,  -1  (40) 

for  2=1  (spherical  coordinates)  or  m=l  (cylindrical  coordinates) . 

For  all  mode  numbers  greater  than  unity  6  is  complex,  and  there  is 

a  factor  growing  with  a  power-low  dependence  ~  r^  (X+a)/(2X)  an(j  a 

1  2 

factor  which  is  oscillatory  in  Jin  r^  with  a  frequency  w  =  (4XQ  “  (^+°0  1 
The  real  part  is  always  negative  and  independent  of  mode  number,  while 
the  oscillatory  part  depends  on  the  mode  number: 


X+q  .  X+a 
C/r o  ^rQ  "  2  X  -  1  X 

where 

[4XQ  -  (X-Kx)  2 1^ 
p  '  2  (X-hx) 


(41) 


(42) 


In  Fig.  2  we  show  the  frequency  in  Jin  t  for  the  spherical  harmonic 

perturbation  as  a  function  of  mode  number  for  two  different  ratios  of 

specific  heat,  V  =  7/5  and  y  =  5/3.  Figure  3  shows  the  ratio  r/r  between 

0 

successive  minima  in  the  perturbation  amplitude  as  it  oscillates  during 
collapse.  The  ratio  may  be  chosen  so  that  at  some  prescribed  degree  of 
compression  the  perturbation  amplitude  is  a  minimum. 


d  A0N3nO3Hd 
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Fig.  2  —  Oscillation  frequency  in  2n  t  of  spherical  harmonic  perturbation  of  a  spherical 
function  of  mode  number  for  a  perfect  gas  with  7=5/3  and  7=7/5. 


COMPRESSION  RATIO  Rn/R 


Fig.  3  —  Compression  ratio  between  successive  minima  in  perturbation  amplitude  during 
perturbed  shock  collapse  as  a  function  of  mode  number. 
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IV.  Numerical  Integration 


In  order  to  assess  the  importance  of  nonlinear  effects  on  the  mode 
amplitude,  the  full  nonlinear  model  equations  were  integrated  with  a  code 
similar  to  that  of  Fong  and  Ahlborn  (1979).  The  code  advances  the  equations 
of  motion  of  the  shock  front 


r  =  v  ; 

dM  _  M2  -  1 
dA  X(M)  M  A  ’ 


(43) 

(44) 


either  in  plane  coordinates  (for  cylindrical  collapse  viewed  in  a  plane 
through  the  axis)  or  in  cylindrical  coordinates  (for  spherical  collapse 
viewed  in  the  equatorial  plane) . 

The  equation  for  the  cross-sectional  area  is  not  integrated  directly, 
but  areas  are  taken  from  the  kinematics  of  the  integrated  ray  tubes.  A 
second-order-accurate  space-and-time-centered  algorithm  is  employed  to 
advance  the  grid  locations  and  Mach  numbers.  Thus  the  problem  of  computing 
the  shock  shape  is  reduced  to  integrating  a  set  of  ordinary  differential 
equations  for  the  trajectories  of  a  finite  number  of  grid  points  located 
along  the  shock  front,  and  for  their  associated  time-dependent  Mach  numbers. 
The  integration  is  subject  to  a  timestep  limit  analogous  to  the  Courant 
condition  for  the  one-dimensional  fluid  equations.  This  algorithm  allows 
for  the  propagation  of  so-called  shock-shocks  (discontinuities  in  the  slope 
and  Mach  number  of  the  shocks  predicted  by  the  Whitham  theory)  in  either 
direction,  since  the  scheme  is  centered  and  symmetric.  The  equations 
for  the  shock  surface  are  analogous  to  the  one-dimensional  Lagrangian 
equation  of  motion,  where  the  shock  position  takes  the  place  of  the  fluid 
position,  the  ray  tube  area  takes  the  place  of  fluid  density,  and  the  Mach 
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number  takes  the  place  of  the  fluid  velocity.  Thus  many  of  the  properties 
of  one-dimensional  motion,  e.g.,  disturbances  traveling  along  characteristic 


directions  and  nonlinear  wave  front  steepening,  show  up  in  the  shape 
of  the  shock  surface. 

It  was  found  necessary  to  redistribute  the  grid  points  to  prevent 
unacceptably  short  timesteps  as  the  mesh  points  crowd  together  near  the 
shock-shock  regions.  The  algorithm  diffuses  the  grid  points  a  small 
amount  parallel  to  the  shock  front  in  such  a  way  as  to  make  the  distances 
between  points  more  nearly  equal.  This  diffusion  plays  a  role  similar 
to  that  of  an  artificial  surface  tension  at  interfaces. 

To  test  the  accuracy  of  the  numerical  procedure  described  above,  we 
compared  the  results  of  the  calculation  with  those  predicted  by  the  self¬ 
similar  solutions  for  the  collapse  of  an  infinitely  strong  shock  due  to 
Guderley  (1942).  The  self-similar  solution  predicts  a  shock  position  given 
ly  R  *  C  (-t)’7  with  V  *  .717  for  ay*  1.4  gas  according  to  Guderley. 


This  gives  a  shock  Mach  number  as  a  function  of  radius 

R  *  -j?  C  (-t)  n_1  ^  R(r'*1)/rl 


Thus 


(46) 


The  C-C-W  approximation  for  large  M  gives 


with 


This  yields 


dA  „  dR  ,,  .  2  l-o2  ,  (2c+l)  M2  +  1  dM 

~  * 2  t  *  (1y+r  ■>!?■- 1 - M 


M  ^  R  (20+1)  (ya+1)  _  -.3941 

5  R 


(47) 


(48) 
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The  exponents  agree  with  the  analytic  Guderley  solution  to  within  0.0015. 
In  Fig.  4  we  show  the  results  of  these  models  for  the  spherical  convergence 
◦f  a  shock  with  an  initial  Mach  number  of  57  at  a  radius  of  0.112  cm. 

The  solid  line  is  the  analytic  result  of  the  Guderley  solution.  The 
triangles  are  the  numerical  results  of  the  c-C-W  approximation  using  the 
code  described  above  with  25  mesh  points  around  a  circle  In  cylindrical 
coordinates  (describing  a  spherical  implosion).  Comparing  the  results 
along  the  axis  with  those  perpendicular  to  the  axis,  we  see  that  the  shock 
remains  spherical  to  better  than  1%  during  the  entire  implosion  and  the 
Mach  number  reproduces  the  Guderley  solution  with  less  than  0.5%  error. 
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Shock  Mach  Number 


Spherical  Shock  Collapse 


Fig.  4  —  Comparison  of  analytic  and  numerical  integration  for  self-similar  spherical  shock 

collapse  (Guderley  problem). 


V.  Comparison  of  Results  of  Linearized  Model  and  Nonlinear  Integration 

In  Fig.  5  and  6  we  show  the  results  of  the  analytic  and  numerical 
integration  for  mode  numbers  2,  4,  and  8  for  cylindrical  and  spherical 
collapse.  The  solid  lines  are  the  analytic  formulas  and  the  circles  show 
the  numerical  results.  To  Insure  beginning  in  the  linear  regime  the  initial 

_3 

amplitude  of  the  modes  was  chosen  such  that  c/r^  “  10  .  Since  the  system 

of  equations  is  second  order  in  the  perturbed  shock  location,  both  the 
amplitude  and  the  velocity  of  perturbation  must  be  specified.  For  these 
cases  the  velocity  was  chosen  to  be  zero.  This  then  uniquely  determines 
the  phase  of  the  oscillation.  The  phase  angle  between  the  amplitude  and 
the  velocity  is  determined  by  the  mode  number  and  the  specific  heat  ratios: 

*-tan  tTTxS)?!  •  (49) 

Agreement  appears  very  good  until  the  mode  amplitude  becomes  greater 
than  a  few  percent.  At  this  time  nonlinear  effects  which  can  generate  modes 
other  than  primary  one  drive  the  solution  away  from  the  linear  result. 

Nonlinear  steepening,  of  the  wave  front  generates  higher  order  modes  which  begin 
to  dominate  the  solution.  In  Fig.  7  we  can  see  the  form  this  takes.  This 
figure  shows  the  successive  wave  front  shapes  for  an  l  ■  8  spherical  har¬ 
monic  perturbation  initialized  with  a  large  initial  amplitude  in  order  to 
show  the  effects  of  a  large  compression.  As  the  shock  collapses,  the 
wave  fronts  begin  to  form  cusplike-shapes  in  the  region  where  the  shock 
is  left  behind.  For  sufficiently  large  amplitudes  a  true  cusp  forms  and 
the  simple  shock  front  no  longer  exists.  It  is  replaced  by  a  system  of 
reflecting  shock  waves  (see  Fig.  8),  which,  however,  are  not  treated  in 
the  present  model.  When  this  occurs  a  large  portion  of  the  shock  energy 
can  be  left  behind  in  the  reflected  shock  system,  thus  decreasing  the 
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compressional  effect  of  the  imploding  shocks.  Solutions  of  this  type  are 
referred  to  below  as  having  kink  instabilities  in  the  shock  front. 
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Fig.  5  —  Comparison  of  analytic  (solid  lines)  and  numerical 
integration  (open  circles)  for  maximum  perturbation  ampli¬ 
tude  as  a  function  of  radius  for  cylindrical  shock  collapse 
at  three  different  mode  numbers.  Deviation  at  small  radii 
is  due  to  nonlinear  effects  not  in  analytic  model. 
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Fig.  6  —  Comparison  of  analytic  (solid  lines)  and  numerical 
results  (open  circles)  for  spherical  shock  collapse  for  mode 
numers  2,  4,  and  8. 
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7  —  An  example  showing  formation  of  cusp-line  structure  for  large  perturbation  ampli- 
3  in  the  nonlinear  regime.  Initial  perturbation  is  a  mode-number-8  Legendre  polynomial 
amplitude  10 %  of  mean  radius. 


Fig.  8  In  the  nonlinear  regime  shock-shocks  are  formed  which  are  the  intersections  of  regular  of  Mach 
reflections  of  shocks.  This  remits  in  potential  loss  of  shock  energy  in  the  reflected  shock  system. 


VI .  Nonlinear  Modeling 


In  order  to  investigate  more  carefully  how  these  kinks  form  and 

propagate,  we  look  at  perturbation  in  the  form  of  a  simple  spherical 

cap  of  radius  smaller  than  the  initial  mean  shock  radius.  This  shock  cap 

intersects  the  main  shock  with  an  angle  6  at  a  colatitude  0,  as  shown 

in  Fig.  9.  For  the  purpose  of  this  investigation  the  Mach  number  around 

the  shock  is  assumed  uniform  at  M  =  10.  In  Fig.  10  we  show  the  limits  of 

stability  for  5  as  a  function  of  8,  In  addition  we  show  the  stability 

result  in  terms  of  Sr^/r^  as  a  function  of  3  (Sr^/r^  is  geometrically 

related  to  6  and  6).  When  either  6  or  0  becomes  too  large,  the  shock  front 

is  transformed  into  a  nonlinear  cusp-type  regime  instead  of  reverting  to 

a  more  spherical  form  and  oscillating.  From  Fig.  10b  we  can  see  that  this 

is  related  to  the  initial  perturbation  and  is  a  function  of  0.  In  each 

case  instability  is  defined  to  occur  if  a  cusp  shock  is  formed  before 
-2 

r/r^<10  .  This  arbitrary  definition  is  necessitated  by  the  fact  that  all 

perturbations,  given  sufficient  compression,  eventually  evolve  into  the 
cusp-shaped  form.  Realistically  speaking,  however,  compressions  of  10^  are 
more  than  sufficient  to  achieve  the  compression  and  temperature  rise 
necessary  for  ignition  in  pellet  fusion. 

Even  for  these  single  cap  perturbations,  oscillatory  behavior  is  apparent 
for  sufficiently  small  initial  amplitudes.  Let  us  compare  the  radius  of 
the  first  minimum  of  the  average  deviation  from  a  spherical  or  cylindrical 
shock  for  the  cap  perturbation  to  that  for  Legendre  polynominal  whose  first 
zero  forms  the  same  angle  as  the  0  for  the  cap  perturbation.  We  see  from 
Fig.  11  that  at  least  for  smaller  0  (i.e.,  larger  l)  the  oscillation  period 
is  nearly  the  same.  This  indicates  that  the  oscillation  period  (i.e.,  the 
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time  or  radius  between  minimum  deviations  from  symmetric  implosion)  can 
be  reasonably  well  predicted  for  arbitrary  perturbations  by  matching  to 
the  lowest  order  Legendre  polynominal  which  fits  the  perturbtation.  In 
Fig.  12  we  show  a  stable  case  where  the  behavior  about  a  nearly  spherical 
implosion  is  oscillatory  up  to  a  compression  of  r^/r  =10.  In  Fig.  13  we 
show  an  unstable  case  where  a  cusp  is  clearly  forming  after  a  single  over¬ 
shoot  of  the  perturbation  and  a  compression  of  r^/r  =  4. 

From  these  results  we  see  that  two  interdependent  factors  control 
when  the  nonlinear  kinks  will  begin  to  form.  One  is  the  wavelength  of 
the  perturbation  and  the  other  is  the  angle  at  which  the  perturbation 
intersects  the  mean  radius.  The  smaller  the  wavelength,  the  larger  the  angle 
that  is  tolerable.  Note,  however,  that  this  angle  and  the  mode  amplitude 
are  not  independent.  In  terms  of  amplitude,  the  shorter  the  wavelength, 
the  smaller  the  amplitude  of  disturbance  that  can  be  tolerated. 
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Fig.  11  —  Compression  ratio  at  the  first  minimum  of  6R/R  as  a  function  of  intersection 
angle.  Open  circles  represent  results  obtained  with  Legendre  polynomials,  where  0  is  the 
angle  to  the  first  zero. 
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Fig.  13  —  Example  of  an  unstable  spherical  cap  perturbation  where  kink  instability  forms  after 
one  oscillation,  (a =39°,  0=34°,  SR/R=0.25) 


VII.  Effect  of  density  variation 

Equation  (I)  was  derived  assuming  a  uniform  density  and  temperature 
ahead  of  the  shocks .  The  C-C-W  methodology  does  not  require  the  medium 
ahead  of  the  shocks  to  be  uniform,  and  we  shall  now  investigate  the  effect 
a  non-zero  density  variation  in  the  undisturbed  medium.  However,  we  do 
not  allow  the  density  to  increase  rapidly  enough  that  the  shock  decelerates. 
This  would  violate  the  important  assumption  that  returning  characteristics 
are  negligible. 

The  compatability  relation  along  the  inward-directed  characteristic 
is  given  by 

2 

dp  +  padu  +  M  =  o  ,  (50) 

while  the  Hugoniot  relations  behind  a  strong  shock  are  given  by 


U  "  Y+l  V 


2  «  2 

p"y+T  pv 

p  ,  XU  p 

p  Y+l  M 


Puv 


2 

a 


JCR 

P 


2V(Y-1) 

(Y+l) 2 


(51) 


where  u  is  the  fluid  velocity  behind  the  shock,  p  is  the  pressure  behind 
the  shock,  p  is  the  density  behind  the  shock,  and  p  is  the  undisturbed 
density.  Differentiating  the  relation  for  p  we  get 


dp 


2 

Y+l 


v2  dp 


p  vdv 


(52) 
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Combining  Eqs.  (50-52)  we  have 


dr  ^  +  <4t  +  VS>  5vdv  +  — 


-  2  dA  (53) 


[2(Y+1)  +  ✓2(Y-1>]  Pv  ~K, 


or 


dv 

v 


-  <*r  +vs>'1  t  - « r 


(54) 


If  we  now  write  this  in  the  form  of  Eq.  (7)  we  have 


X  (v/v)  ”  -A/A  -  <  p/p  . 


(55) 


Using  Eqs.  (4),  (6)  and  (10)  and  the  relation 


P  ■  v  •  V  p  , 


we  have 


(56) 


(X/v)  *V*n+<n»7JlnP 


(57) 


Assuming  a  density  profile  of  the  form  p  -  r~q  the  zeroth  order 
symmetric  equation  becomes  similar  to  Eq.  (15), 


(X/r)*-  (a-K)/r 


(58) 
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This  has  a  solution  analogous  to  Fq.  (18), 


r 


R 


(59) 


The  effect  of  a  positive  q  (density  increasing  toward  the  axis)  is  to  cause 
the  shock  to  accelerate  more  slowly  than  the  q“0  case.  (A  negative  q 
will  cause  the  shock  to  accelerate  faster). 

Again  linearizing  the  equations  as  in  Section  2,  we  have  to  first 

order 


-  V-i "  (7oP :  ( W  +  f  <  ni‘  Vn  5 

-"o  *(V  ‘  70  zn  P  +  "o  *  70  *  70  ln  ^ 


(60) 


Using  vector  manipulations  similar  to  those  in  Section  2,  we  obtain 


-(XnQ*C  /v“)  -  -  ?Q*  [(I-n0n0)»  V  (nQ*S))  +  ^  *  v0  (®/ro)  + 


(61) 


-o2  :  Vo  ln  p 


The  last  term  of  Eq.  (61)  adds  an  additional  term  to  Eq.  (29),  given  by 


-0  VO  ln  p  m  C(£np)  * 


(62) 


where  the  double  prime  indicates  a  derivative  with  respect  to  r^. 
Thus  Eq.  (61)  becomes 
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(63) 


We  note  from  this  result  that  the  effect  of  including  density  variation 
is  found  to  first  order  by  replacing  everywhere  the  term  a  by  m-fcq. 

That  is,  the  acceleration  normally  due  to  the  geometric  factor  <x  is  modified 
by  an  additional  term  < q  related  to  the  acceleration  caused  by  a  density 
variation.  Positive  q  (density  increasing  toward  the  axis)  has  the  effect 
of  reducing  the  acceleration  due  to  geometric  convergence.  For  sufficiently 
large  q  the  acceleration  vanishes.  This  situation,  however,  violates  the 


assumption  in  the  C-C-W  model  that  the  shock  system  is  accelerating.  The 
ratio  of  the  disturbance  amplitude  to  the  2eroth-order  radius  continues  to 


satisfy  Q/r 


0 


,  with  the  frequency  appropriately  modified.  In  terms 


of  the  zeroth  order  radius,  however,  the  growth  becomes  C/r  ~  r, 


-(X+a-tcq)  /  (2X) 


‘0  ‘0 

A  positive  q  thus  decreases  the  relative  growth  of  the  shock  perturbation 
as  a  function  of  the  shock  radius. 
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VHI.  Conclusions 


The  C-C-W  approximation,  which  appears  to  be  very  accurate  for  com-  II 

puting  converging  shock  waves,  predicts  that  the  converging  shock  is  always  9 

unstable,  in  the  sense  that  the  ratio  of  the  perturbation  size  to  the  H 

average  radius  diverges  at  the  time  of  collapse.  The  cylindrical  and  9 

spherical  cases  differ  only  quantitatively.  The  growth  rate  is  only  a  9 

power  law,  however,  and  therefore  is  not  as  serious  as  an  exponential  growth.  X| 

For  mode  number  greater  than  unity,  the  perturbations  oscillate  in  In  t 
with  a  mode-number-dependent  period.  The  amplitude  growth,  however,  is 
independent  of  mode  number. 

For  sufficiently  large  amplitude  the  linear  behavior  breaks  down 
and  the  solution  develops  nonlinearly  into  a  kink  form,  where  a  reflected 
shock  is  left  behind  in  the  material .serving  as  a  potential  loss  of  shock 
energy  in  the  collapsing  shock.  [ 
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